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Abstract. We present results from numerical simulations of Ray leigh- Taylor 
turbulence, performed using a recently proposed [TJ [2] lattice Boltzmann method able 
to describe consistently a thermal compressible flow subject to an external forcing. 
The method allowed us to study the system both in the nearly-Boussincsq and strongly 
compressible regimes. Moreover, we show that when the stratification is important, 
the presence of the adiabatic gradient causes the arrest of the mixing process. 
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1. Introduction 

Computational methods based on discrete-velocity models have gained considerable 
interest in the recent period as efficient tools for the theoretical investigation of the 
properties of complex flows [3J [4J, [5J [6j [T] . In particular, it has been recently shown that an 
important class of these models, known as the lattice Boltzmann models (LBM) [8|l9llT0]. 
can be derived from the continuum Boltzmann (BGK) equation [11]. This derivation 
involves the expansion in suitable Hermite polynomials of the distribution functions 
f(x,£,t), describing the probability of finding a molecule at space-time location (x,t) 
and with velocity £ [51 [121 E2 [II]- Therefore, the corresponding lattice dynamics is 
well founded in terms of an underlying continuum kinetic theory. The state-of-the-art 
is satisfactory concerning iso-thermal flows, even in presence of complex bulk physics 
(multi-phase, multi-components) [31 HI [15] and/or with complex boundary properties 
including roughness, wetting and slip boundary conditions [161 EL71 H] • 
The situation is much less satisfactory when hydrodynamical temperature fluctuations 
play an active role in the flow evolution, due to complex compressible effects or to 
phase-transition in multi-phase systems. 

Within this framework, we have recently developed [II [2] a new LBM which allows to 
incorporate the effects of external/internal forces in thermal systems. Here we use this 
new algorithm to study highly compressible Rayleigh- Taylor systems, with an initial 
configuration such that two blobs of the same fluid are prepared with two different 
temperatures (hot, less dense, blob below, cold, denser, blob above). We show that the 
method is able to handle the highly non-trivial spatio-temporal evolution of the system 
even in the developing turbulent phase. In this case, we could push the numerics up to 
Atwood numbers At ~ 0.4. Maximum Rayleigh numbers achieved are Ra ~ 4 x 10 10 for 
At = 0.05 and Ra ~ 2 x 10 9 for At = 0.4. The paper is organized as follows. We will first 
describe the method (section [2]), the numerical setup (section [3]) and the system studied 
(section H|); then we will discuss the two main physical results, i.e. the stratification 
(section |5]) and compressibility (section [6]) effects, and some features related to the 
conservation of mean quantities (section [7]). Conclusions are in section [HJ 

2. The Lattice Boltzmann Model 

Here we introduce the simulated equations set along with a brief description of the 
computational lattice Boltzmann method employed. More details, along with many 
validations, can be found in [2]. The thermal-kinetic description of a compressible 
gas/fluid of variable density p, local velocity u, internal energy /C, and subject to a 
local body force density g, is given by the following equations: 

d t p + di(pUi) = 
d t (pu k ) + di{P ik ) = pg k 

d t fC + -dM = pgiUi (1) 
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Figure 1. Scheme of the discrete set of velocities; r is the lattice constant whose value 
is r = 1.1969... (see [2] and references therein). 



where Pik and qi are the momentum and energy fluxes, describing advection, viscous 
properties and thermal diffusivities in the hydrodynamical limit. 

In pQ it has been shown that it is possible to recover exactly these equations, 
starting from a continuum Boltzmann Equation and introducing a suitable shift of the 
velocity and temperature fields entering in the local equilibrium. 

The lattice counterpart of the continuum description can be obtained through the 
lattice Boltzmann discretization (fi(x,t) are the fields associated to the populations): 

fi(x + Cl At,t + At) - /,(*,*) = ~ (/,(*,*) - fl eq) ) (2) 

where the equilibrium is expressed in terms of hydrodynamical fields and the body 
force term g, and the subscript / runs over the discrete set of velocities q (see fig. [TJ; 
in equation (|2J) r is the relaxation time (which is related to the dynamic viscosity v via 
v = pT(r — 1/2), T being the temperature field), and At the time step of the simulation. 

The macroscopic fields (density, momentum and temperature) are defined in terms 
of the lattice Boltzmann populations: p = Y^i fi, pu = Y^i cifi, DpT = ^ \c t — u\ 2 //. 
Lattice discretization induces non trivial corrections terms in the macroscopic evolution 
of averaged hydrodynamical quantities. In particular, both momentum and temperature 
must be renormalized by discretization effects in order to recover the correct description 
out of the discretized LBM variables: the first correction to momentum is given by 
the pre and post-collisional average PU ED] = u + Ai g an d t he first, non-trivial, 
correction to the temperature field by [T] = T + ■ (D is the dimensionality of 

the system). Using these "renormalized" hydrodynamical fields it is possible to recover, 
through a Taylor expansions in At, the thermo-hydrodynamical equations [H[2]: 

D t p = -pd^ (3) 
pD t ul H) = -dip - pg5 i)3 + udjjU^ (4) 
pc v D t T^ + pd t uf ] = kd u TW (5) 

where we have introduced the material derivative, D t = dt + U- dj, and we have 
neglected viscous dissipation in the temperature equation (usually small). Moreover, 
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Table 1. Parameters for the three types of Rayleigh- Taylor runs. Atwood number, 
At = (Td —T u )/ (Td + T U ); viscosity v; gravity g; temperature in the upper half region, 
T u ; temperature in the lower half region, T^; normalization time, f = yj L x /(g At). 

c v is the specific heat at constant volume for an ideal gas p = pT^ H \ and v and k are 
the transport coefficients. From now on, for the sake of simplicity, we will drop the 
superscript (H), knowing that we are dealing with lattice hydrodynamical quantities 
satisfying equations (J3H5]). As a tool for the numerical simulation of systems as (or 
similar to) the one we plan to study, the LBM may suffer, in principle, of some issues, 
such as having too high Mach numbers, too low viscosity (i.e. very small relaxation 
time, which is undesirable especially in the presence of processes very far from local 
equilibrium): we could, however, check the accuracy of our method against different 
ones, finding it extremely competitive, within the range of parameters discussed in this 
paper [2~T] . 

3. Details of the numerical simulations 

We use a 2D LBM algorithm, with 37 population fields (a so called D2Q37 model), 
moving with the lattice speeds shown in figure (CQ). Since the lattice spacing can be 
taken to be unitary, the time step At will be the inverse of the lattice unit speed, 
i.e. At ~ 0.835.... Three different kinds of simulations have been performed (whose 
parameters are summarized in table CD): (A) with a large enough adiabatic gradient 
(but small Atwood number) in order to address the stratification effects on the mixing 
layer growth, while still being very close to the Boussinesq approximation; (B) with an 
adiabatic gradient which is twice the one of run A; (C) with large Atwood in order to 
describe compressibility effects, out of the Boussinesq regime, but far from the adiabatic 
profile; (D) with small adiabatic gradient and small Atwood number. 

4. The Rayleigh- Taylor system 

Superposition of a heavy fluid above a lighter one in a constant acceleration field depicts 
a hydrodynamic unstable configuration called the Rayleigh- Taylor (RT) instability 
[22] with applications on different fields going from inertial-confinement fusion [23] 
to supernovae explosions [21] and many others [25]. Although this instability was 
studied for decades it is still an open problem in several aspects [26]. In particular, 
it is crucial to control the initial and late evolution of the mixing layer between the two 
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Figure 2. Initial condition of the Rayleigh- Taylor system, as given by equation Q. 
We show the mean temperature (bold line) and density (tiny line) profiles as function 
of the z-coordinate (on the vertical axis). 



miscible fluids; the small-scale turbulent fluctuations, their anisotropic/isotropic ratio; 
their dependency on the initial perturbation spectrum or on the physical dimensions of 
the embedding space [27J EE] • In many cases, especially concerning astrophysical and 
nuclear applications, the two fluids evolve with strong compressible and/or stratification 
effects, a situation which is difficult to investigate either theoretically or numerically. 
Here, we concentrate on the large scale properties of the mixing layer, studying a slightly 
different RT system than what usually found in the literature: the spatio temporal 
evolution of a single component fluid when initially prepared on the hydrostatic unstable 
equilibrium, i.e. with a cold uniform region in the top half and a hot uniform region 
on the bottom half (see figure [2]). For the sake of simplicity we limit the investigation 
to the 2d case. While small-scales fluctuations may be strongly different in 2d or 3d 
geometries, the large scale mixing layer growth is not supposed to change its qualitative 
evolution (29J [30]. Grey-scale coded snapshots of a typical RT evolution are shown in 
figure |3] showing all the complexity of the phenomena. Let us start to define precisely 
the initial set-up. We prepare a single component compressible flow in a 2d tank of size, 
L x x L z , with adiabatic and no-slip boundary conditions on the top and bottom walls, 
and with periodic boundary conditions on the vertical boundaries. For convenience we 
define the initial interface to be at height z = 0, the box extending up to z = L z /2 
above and z = —L z /2 below it (see figure |2]). In the two half volumes we then fix two 
different homogeneous temperature, with the corresponding hydrostatic density profiles, 
p , verifying 



Considering that in each half we have po(z) = Tpo(z), with T fixed, the solution has an 
exponentially decaying behavior in the two half volumes, each one driven by its own 
temperature value. The initial hydrostatic unstable configuration is therefore given by: 




(6) 



T (z) = T u ; p (z) = p u exp(-g(z - z c )/T u ); 



z>0 
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Figure 3. Spatio-temporal evolution for a typical RT run with L x x L z = 1024 x 2400, 
T u = 0.95, Td = 1.05 at four instant of time: t = f, 2f , 4f , 6f (run D in table Q} going 
clockwise from the top left panel. 



ToO) = T d \ p (z) = p b exp(-g(z - z c )/T d ); z < 0. 

(7) 

To be at equilibrium, we require to have the same pressure at the interface, z = z c = 0; 
which translates in a simple condition on the prefactor of the above expressions: 

PuT u = PbT d . (8) 

Because T u < T d , we have at the interface p u > p b . As far as we know, there are 
no exhaustive detailed calculations of the stability problem for such configuration, even 
though not too different from the usual RT compressible case [22} [HI [32] . As said, this is 
not the common way to study RT systems, which is usually meant as the superposition 
of two different miscible fluids, isothermal, with different densities [22j [331 EB I2Z]- As 
far as compressible effects are small, one may safely neglect pressure fluctuations and 
write - for the case of an ideal gas: 
5p 5T 
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and the two RT experiments are then strictly equivalent. Moreover, in the latter case, if 
one may neglect the dependency of viscosity and thermal diffusivity from temperature, 
the final evolution is undistinguishable from the evolution of the temperature in the 
Boussinesq approximation [29| [28] . 



5. The adiabatic gradient and the arrest of the mixing process 

The main novelty in the system here investigated is due to the presence of new effects 
induced by the adiabatic gradient, which can be written for an ideal gas as (3 a d = g/c p . 
The role of stratification, i.e. of the adiabatic gradient, is quite well established in 
the context of Rayleigh- Benard convection (see, e.g. [S]), while it has been only in 
recent times studied, numerically [351 EH] and theoretically [37], in a set up such as 
that of Rayleigh- Taylor mixing. In order to understand the main physical point it is 
useful to think at the RT mixing layer as equivalent to a (developing) Rayleigh-Benard 
system with an imposed mean temperature gradient [381 [39]. Let us denote with L m i(t) 
the typical width of the RT mixing layer at a given time as measured for example 
from the distance between the two elevations where the mean temperature profile is 
1% lower or higher then the bottom and top, respectively, unmixed temperature values, 
L m i = z u — z d , where (T(x,z u )) x = 1.01T U and (T(x,Zd)) x — 0.99T d . The temperature 
tends to develop a linear profile inside the mixing region, the resulting instantaneous 
temperature gradient being given by j3{t) = {T d — T u ) / L m i(t) , and, hence, decreasing in 
time inversely to the growth of the mixing length. As a result, at a certain time (if the 
box is high enough) the instantaneous temperature gradient will become of the same 
order of the adiabatic gradient, (3(t) ~ (3 a d and the growth of the mixing length will 
stop. In figure (j4j) we show the mean temperature profiles once the mixing is already 
stopped, for two different values of gravity (runs A and B in table (CQ)): in the mixing 
layer the two curves have developed a linear profile with slope g/c p , which is exactly the 
adiabatic gradient for an ideal gas. One can define an instantaneous Rayleigh number, 
driving the physics inside the mixing layer, as: 

~ (g/f )LUt)(m-Pa d ) nm 
na{t) = — — — , (1U) 

(k/p c p ){v/p ) 

where (•) indicates quantities evaluated at the middle layer. It is clear that for small 
times, (3(t) <C (3 ad , the effective instantaneous Rayleigh number is high: the system is 
unstable, and the mixing length grows. On the other hand, as time elapses, the vertical 
mean temperature gradient decreases, until a point when, (3(t) ~ f3 ad , the instantaneous 
effective Rayleigh number becomes Ra(t) ~ 0(1) and the system tends to be stabilized. 
We can then identify an adiabatic length: 

Lad = (Td - T u )/f3 ad = c p AT/g 



which determines the maximum length achievable by the mixing layer, in our 
configuration. When the mean temperature approaches the adiabatic profile, the system 
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Figure 4. Mean temperature profiles (T) x (z,t) for run A (A) and run B (0), in 
a stage where the mixing process is already stopped (t ~ 13f for both cases). The 
dashed lines represent the correspoding adiabatic profiles. 



shows a sudden slowing down of the mixing layer growth which, eventually, stops. A 
possible way to estimate quantitatively when and how the adiabatic gradient starts to 
play a role in the growth of the mixing length is to use a simple phenomenological 
closure for large scale quantities in the system. We start from the self-similar scaling 
predicted by [10, HI] for the homogeneous unstratified growth: 

(L ml (t)) 2 = Aa g At L ml (t) (11) 

which has a unique solution in terms of the initial value, L m i(to): 

L ml {t) = L ml (t ) + 2y/aAtg(t-t ) + aAtg(t- t ) 2 . (12) 

In order to minimally modify the above expression considering the saturation effects 
induced by stratification, we proposed to use in [2]: 

(L ml (t)) 2 = AagAtL^m (h^l) (13) 



where ip = ip{x) must be a function fulfilling the condition ip — y 1 as x — > (that is 
for L ad — y oo), in order to recover the equation ffTTj) for the unstratified case when the 
adiabatic gradient goes to zero. We further add the requirement of reaching the adiabatic 
profile with zero velocity and acceleration, enforcing a strict irreversible growth, i.e. 
L m i > 0, as it must be for the case of miscible fluids. Under these assumptions it can 
be shown that the simplest form for the function if) is: 

where the prefactor C must be set equal to l/(e — 2) to comply with the 
prescribed boundary conditions. Equation ( 113]) must be considered as a zero-th order 
phenomenological way to take into account for the adiabatic gradient in the mixing 
layer evolution. We integrated numerically eq. ( TT3]) testing the result in figure [5] where 
we show that it is possible to fit the global evolution of the mixing length L m i(t), by 
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Figure 5. Evolutions of the mixing layer, L m i(t) versus time with two different 
adiabatic lengths, corresponding to run A (A) and run B (0) in table [TJ Solid bold 
lines correspond to the theoretical prediction (fT3f with a = 0.05. 



using reasonable [26] values of a, for all times, including the long time behavior where 
Lmiif) ~ L a d. We can then interpret the solution of our equation ( Tl3l) as a good 
generalization of f fl2|) including also the adiabatic gradient effects. 

6. Effects of compressibility 

In this section we are going to study the effects of flow compressibility on the dynamics 
of a RT system at changing Atwood number. To do that we come back to the discussion 
sketched in section HJ since the equation of state for our fluid is the one of a perfect gas, 
the pressure, temperature and density fluctuations are related by: 
5P 5p 5T 

T = j + Y' (15) 

hence, as discussed in section HJ if pressure fluctuations are small, density fluctuations 
are linearly slaved on those of temperature, which will be also small, and the system 
behaves as a Boussinesq fluid. Conversely, if the temperature jump is high there will be 
large density differences through the system ([7]) and hence large pressure fluctuations 
((6]); thus we expect that, at increasing the Atwood number, the dynamics becomes more 
and more compressible and pressure turns out to be a dynamically relevant variable. 
We show, first, how the mixing acts on the statistics of the (point-wise) temperature; as 
we can see in figures ([6]) and ([7]), where we plot the probability density function (PDF) 
of temperature at two instants of time for the two At = 0.05 and At = 0.4 (for which the 
squared Mach number is Ma 2 ~ 0.16), the distribution has initially a bimodal character, 
since at the beginning the volume is divided into two homogeneous regions of hot and 
cold fluid. Due to the mixing, at later times, the probability of having intermediate 
values of temperature increases; however, the two peaks remain dominant, because the 
system dynamics does not involve yet the whole box and because the diffusive processes 
are so slow to be irrelevant at this stage. No evident differences (except the obviously 
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Figure 6. Probability density functions of temperature at t = f (©) and at t = 5f 
(A), for At = 0.05, for which the flow is basically Boussinesq-likc. 
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Figure 7. Probability density functions of temperature at t = f (©) and at t = 5f 
(A), for At = 0.4. that is in the highly compressible regime. 

larger range of values spanned) rise between the low and high Atwood number cases. 
Hence, to better address this point, we study the statistics of pressure fluctuations (with 
respect to the mean profile), in the two compressibility regimes; we define the fluctuation 
of the generic thermohydrodynamic field <f> as 

<ty(x,t) = 0(x,t) - {<j>) x (z,t)\ ((j)) x (z,t)= (f)dx. 

Jo 

In figure ([H]) we show the PDFs of Sp, measured again inside the whole volume, at two 
different instants of time during the mixing process, for two Atwood numbers. It can 
be noticed how the PDF, while being basically a 5-function for low At and remain- 
ing such at any time (see inset of figure ©), is more spread at higher At and enlarges 
its tails as the time elapses, confirming that pressure dynamics is now highly non-trivial. 

It is, moreover, known that increasing the degree of compressibility of the dynamics 
has also a strong impact in the stability properties of the system (12] and in the statistics 
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Figure 8. Probability density functions of pressure fluctuations at time t = f (©) and 
t = 5f (A) for At = 0.4 and (inset) for At = 0.05 (t — f, solid line, and t = 5f , dashed 
line): while in the low At case the PDF remains basically a 5— function at any time, it 
is more spread (with tails becoming larger as the time increases) in the compressible 
case (At = 0.4). 



of the mixing layer growth process, determining in the latter case an asymmetry in 
the growth of the mixing layer [32], noticeable also in the statistics of the growth 
parameter a [2]. We would like to discuss here such effects, without appealing to 
any phenomenological model, but in terms, again, of PDFs of the temperature field. 
The use of PDFs to address compressibility effects in RT was also suggested, although 
in a slightly different way, in [J3], in regimes from low (Ma 2 ~ 0.008) to moderately 
high (Ma 2 ~ 0.1) squared Mach number. With this aim, we measured the V(T) where 
T = T(x, z*,t* = 5f) along lines at two fixed heights z* (at a certain time in the late 
stage of the evolution), within the mixing layer, symmetrically displaced with respect 
to the mid cell; in particular we chose z* = ±L m i(t)/2. In figures and (flOl) we plot 
the PDFs, for such heights, for At = 0.05 and At = 0.4 respectively. In bothe cases, of 
course, the PDF corresponding to the upper height shows a peak at lower values of T 
(close to the unmixed cold fluid value), and vice versa for the lower height. However, 
while for small At the two PDFs are symmetric with respect to the average temperature 
(in some sense they transform into each other upon reversal around T = 1), for the 
compressible case figure (flOl) displays a clear asymmetry, where the PDF measured at 
the lower z— location develops a more intense tail at low T values, indicating that falling 
cold fluid spikes are faster (and mix slowly) than rising hot fluid bubbles. 

7. Evolution of averaged quantities 

If we integrate equation @, multiplied bu Ui, over the whole volume, since the boundary 
conditions are periodic at the vertical walls and set zero velocity (no-slip) at top and 
bottom plates, we obtain the following equation for the mean kinetic energy: 

pu 2 

d *(~2~) = (P9 U z) ~ e diss (16) 
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Figure 9. PDFs of the temperature field T(x,z*,t = 5f), measured at z* = 
±L m i(t)/2, for At — 0.05. The two PDFs show peaks close to the values of T of the 
unmixed cold (at z* = +L m i/2, (A)) and hot (at z* = —L m i/2, (©)) fluid respectively; 
the two PDFs are symmetric to each other with respect to the mean temperature value 
T = 1 (typical phenomenology for a Boussinesq system) . 





Figure 10. Same as figure ([9]) but for At = 0.4. Differently from the Boussinesq (low 
At) case, the two PDFs are not symmetric any more, since the one measured at the 
lower height develops a fatter tail at low T values, indicating the asymmetry in the 
mixing process evolution. 



where (...) = l/(L x L z ) J(...)dxdz denotes the space average, and edi SS = v((djUi) 2 ). 
Equation ( TT6|) indicates that the total forcing, due to the gravitational field, is consumed 
partly being transformed in kinetic energy and partly by dissipation. In figure (fTTT) 
we show the fraction of forcing transferred as kinetic energy ((dE/dt)/(pgu z )) and 
dissipated (e^ ss / (pgu z )), and we observe that the latter is much smaller than the former, 
as one would have expected, the dissipation being negligible in 2D. This behaviour is in 
strong contrast with what happens in 3D, where energy is transferred downscales and a 
sort of equipartition between kinetic energy growth and energy dissipation is achieved 
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Figure 11. Time derivative of the kinetic energy (©) and dissipation (A) as functions 
of time, shown as fraction of the forcing term (pgu z ). Notice that the forcing due to 
gravity is almost entirely converted in kinetic energy, since the dissipation is basically 
negligible in 2D turbulence. 



8. Conclusions 

We simulated, by means of a new lattice Boltzmann algorithm, the turbulent dynamics 
of a Rayleigh- Taylor system, the characteristics of the method letting us tune the effects 
of both stratifications and compressibility. Concerning the former problem, we discussed 
the importance of the adiabatic gradient for the growth of the RT mixing layer, showing 
the existence of the phenomenon of the arrest of the mixing process and of a maximal 
width, the adiabatic length, L a d, for the mixing region. We measured, then, the PDFs of 
density and temperature fluctuations, inside the mixing region, observing that while the 
two statistics are almost identical for small Atwood numbers (negligible compressibility), 
they decouple when compressibility is large, owing to the increased relevance of pressure 
fluctuations, whose PDF we also measured, confirming that for large temperature jumps 
pressure plays an active dynamical role. 
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